### Code for "Democratization and Economic Output in Sub-Saharan Africa"
### Daniel De Kadt and Stephen B. Wittels

## A note to users:
## Use setwd() to set the working directory to the location where data files are saved. 
## Figures are programmed to be saved automatically to the working directory.
## They are named according to their figure number in the paper.
## The three tables in the paper are saved as the objects "mali.weights," "panel.estimates," and
## "moderators." Code to print these objects to the console is included at the end.


## Options and Libraries
options(scipen = 6, digits = 3)

library(foreign)
library(Synth)
library(xtable)
library(rgenoud)
library(reshape2)
library(quadprog)
library(ucminf)
library(Rcgmin)
library(Rvmmin)
library(minqa)
library(Rcpp)
library(ggplot2)
library(plyr)
library(grid)
library(lme4)


## Data
load("afripanel_wdk_final.RData")
a <- read.csv("conditioning_variables1.csv")
panel.reg <- read.dta("panel.reg1.dta")


## Functions
synth.wdk1 <- function(
  unitID,
  fullname,
  begin,
  end,
  tr2,
  final,
  low,
  high
){
  
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  controls <- unique(data$WBCode[data$WBCode!=unitID])
  
  prep <- dataprep(
    foo=data, 
    predictors=c(
      "lngdpmadlag",
      "lngdpmadlag2",
      "lngdpmadlag3",
      "lngdpmadlag4",
      "lnpop",
      "ki",
      "openk",
      "civwar",
      "civwarend",
      "pwt_xrate",
      "pwt_xrate_lag1",
      "pwt_xrate_lag2",
      "pwt_xrate_lag3",
      "wbank",
      "wbank_lag1",
      "wbank_lag2",
      "imfadj",
      "imfadj_lag1"
    ),
    dependent="lngdpmad",
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out <- synth(prep)
  
  path.plot(synth.res=out, dataprep.res=prep,
            Ylab="Log GDP per capita", Legend=NA, tr.intake=tr2,
            Ylim=c(low,high) , Main=fullname
  )
  
  gap <- prep$Y1plot - (prep$Y0plot %*% 
                          out$solution.w)
  
  colnames(gap) <- unitID
  
  synth.table <- synth.tab(dataprep.res=prep, synth.res=out)
  
  print(xtable(synth.table$tab.w))
  
  print(xtable(synth.table$tab.v))
  
  tab <- cbind(prep$X1,prep$X0%*%out$solution.w)
  colnames(tab) <- c("Treated Unit","Synthetic Unit")
  print(tab)
  
  return(gap)
}

synth.wdk2 <- function(
  unitID,
  fullname,
  begin,
  end,
  tr2,
  final,
  low,
  high
){
  
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  controls <- unique(data$WBCode[data$WBCode!=unitID&data$WBCode!="ETH"&data$WBCode!="SDN"])
  
  prep <- dataprep(
    foo=data, 
    predictors=c(
      "lngdpmadlag",
      "lngdpmadlag2",
      "lngdpmadlag3",
      "lngdpmadlag4",
      "lnpop",
      "ki",
      "openk",
      "civwar",
      "civwarend",
      "pwt_xrate",
      "pwt_xrate_lag1",
      "pwt_xrate_lag2",
      "pwt_xrate_lag3",
      "eximdiff",
      "eximdiff_lag1",
      "eximdiff_lag2",
      "wbank",
      "wbank_lag1",
      "wbank_lag2",
      "imfadj",
      "imfadj_lag1",
      "imfadj_lag2"
    ),
    dependent="lngdpmad",
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out <- synth(prep)
  
  path.plot(synth.res=out, dataprep.res=prep,
            Ylab="Log GDP per capita", Legend=NA, tr.intake=tr2,
            Ylim=c(low,high) , Main=fullname
  )
  
  gap <- prep$Y1plot - (prep$Y0plot %*% 
                          out$solution.w)
  
  colnames(gap) <- unitID
  
  synth.table <- synth.tab(dataprep.res=prep, synth.res=out)
  
  print(xtable(synth.table$tab.w))
  
  print(xtable(synth.table$tab.v))
  
  tab <- cbind(prep$X1,prep$X0%*%out$solution.w)
  colnames(tab) <- c("Treated Unit","Synthetic Unit")
  print(tab)
  
  return(gap)
}

synth.wdk2.paths <- function(
  unitID,
  fullname,
  begin,
  end,
  tr2,
  final,
  low,
  high
){
  
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  controls <- unique(data$WBCode[data$WBCode!=unitID&data$WBCode!="ETH"&data$WBCode!="SDN"])
  
  prep <- dataprep(
    foo=data, 
    predictors=c(
      "lngdpmadlag",
      "lngdpmadlag2",
      "lngdpmadlag3",
      "lngdpmadlag4",
      "lnpop",
      "ki",
      "openk",
      "civwar",
      "civwarend",
      "pwt_xrate",
      "pwt_xrate_lag1",
      "pwt_xrate_lag2",
      "pwt_xrate_lag3",
      "eximdiff",
      "eximdiff_lag1",
      "eximdiff_lag2",
      "wbank",
      "wbank_lag1",
      "wbank_lag2",
      "imfadj",
      "imfadj_lag1",
      "imfadj_lag2"
    ),
    dependent="lngdpmad",
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out <- synth(prep)
  
  path.plot(synth.res=out, dataprep.res=prep,
            Ylab="Log GDP per capita", Legend=NA, tr.intake=tr2,
            Ylim=c(low,high) , Main=fullname
  )
  
  gap <- prep$Y1plot - (prep$Y0plot %*% 
                          out$solution.w)
  
  colnames(gap) <- unitID
  
  synth.table <- synth.tab(dataprep.res=prep, synth.res=out)
  
  print(xtable(synth.table$tab.w))
  
  print(xtable(synth.table$tab.v))
  
  tab <- cbind(prep$X1,prep$X0%*%out$solution.w)
  colnames(tab) <- c("Treated Unit","Synthetic Unit")
  print(tab)
  
  paths <- cbind(prep$Y1plot,prep$Y0plot%*%out$solution.w)
  
  return(list(paths,out))
}


placebo.model.check <- function(
  unitID,
  begin,
  end,
  tr2,
  final){
  
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  controls1 <- unique(data$WBCode[data$WBCode!=unitID])
  
  prep1 <- dataprep(
    foo=data, 
    predictors=c(
      "lngdpmadlag",
      "lngdpmadlag2",
      "lngdpmadlag3",
      "lngdpmadlag4",
      "lnpop",
      "ki",
      "openk",
      "civwar",
      "civwarend",
      "pwt_xrate",
      "pwt_xrate_lag1",
      "pwt_xrate_lag2",
      "pwt_xrate_lag3",
      "wbank",
      "wbank_lag1",
      "wbank_lag2",
      "imfadj",
      "imfadj_lag1"
    ),
    dependent="lngdpmad",
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls1,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out1 <- capture.output(synth(prep1))
  out1 <- as.numeric(gsub("[A-Z :)()]","",out1[13]))
  
  controls2 <- unique(data$WBCode[data$WBCode!=unitID&data$WBCode!="ETH"&data$WBCode!="SDN"])
  
  prep2 <- dataprep(
    foo=data, 
    predictors=c(
      "lngdpmadlag",
      "lngdpmadlag2",
      "lngdpmadlag3",
      "lngdpmadlag4",
      "lnpop",
      "ki",
      "openk",
      "civwar",
      "civwarend",
      "pwt_xrate",
      "pwt_xrate_lag1",
      "pwt_xrate_lag2",
      "pwt_xrate_lag3",
      "eximdiff",
      "eximdiff_lag1",
      "eximdiff_lag2",
      "wbank",
      "wbank_lag1",
      "wbank_lag2",
      "imfadj",
      "imfadj_lag1"
    ),
    dependent="lngdpmad",
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls2,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out2 <- capture.output(synth(prep2))
  out2 <- as.numeric(gsub("[A-Z :)()]","",out2[13]))
  
  res <- ifelse(out1<out2,1,2)
  
  return(res)
}

synth.placebo1 <- function(
  unitID,
  begin,
  end,
  tr2,
  final
){
  
  controls <- placebo[which(placebo!=unitID)]
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  
  prep <- dataprep(
    foo=data, 
    predictors=c("lngdpmadlag",
                 "lngdpmadlag2",
                 "lngdpmadlag3",
                 "lngdpmadlag4",
                 "lnpop",
                 "ki",
                 "openk",
                 "civwar",
                 "civwarend",
                 "pwt_xrate",
                 "pwt_xrate_lag1",
                 "pwt_xrate_lag2",
                 "pwt_xrate_lag3",
                 "wbank",
                 "wbank_lag1",
                 "wbank_lag2",
                 "imfadj",
                 "imfadj_lag1"
    ),
    dependent="lngdpmad", 
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out <- synth(prep)
  
  gap <- prep$Y1plot - (prep$Y0plot %*% 
                          out$solution.w)
  colnames(gap) <- unitID
  
  return(gap)
}

synth.placebo2 <- function(
  unitID,
  begin,
  end,
  tr2,
  final
){
  
  controls <- placebo[which(placebo!=unitID&placebo!="ETH"&placebo!="SDN")]
  data <- afripanel[afripanel$WBCode==unitID | afripanel$cont_dem_ind==1,]
  
  prep <- dataprep(
    foo=data, 
    predictors=c("lngdpmadlag",
                 "lngdpmadlag2",
                 "lngdpmadlag3",
                 "lngdpmadlag4",
                 "lnpop",
                 "ki",
                 "openk",
                 "civwar",
                 "civwarend",
                 "pwt_xrate",
                 "pwt_xrate_lag1",
                 "pwt_xrate_lag2",
                 "pwt_xrate_lag3",
                 "eximdiff",
                 "eximdiff_lag1",
                 "eximdiff_lag2",
                 "wbank",
                 "wbank_lag1",
                 "wbank_lag2",
                 "imfadj",
                 "imfadj_lag1",
                 "imfadj_lag2"
    ),
    dependent="lngdpmad", 
    unit.variable="wbcode2",
    time.variable="year", 
    treatment.identifier=unitID, 
    controls.identifier=controls,          
    time.predictors.prior=c(begin:end),
    time.optimize.ssr=c(begin:tr2),
    time.plot=c(begin:final),
    unit.names.variable="WBCode"
  )
  
  out <- synth(prep)
  
  gap <- prep$Y1plot - (prep$Y0plot %*% 
                          out$solution.w)
  colnames(gap) <- unitID
  
  return(gap)
}

## Figure 1 -- Left panel
mali <- synth.wdk2("MLI", "", 1975, 1990, 1991, 2008, 6, 8)  
abline(v=1992,lty=3, lwd=2,col="gray")    ## First election
legend("topleft", bty="n", legend=c("Mali", "Synthetic Mali"),
       lty=c(1,2),lwd=c(2,2), cex=1.2)

pdf(file="fig1-left.pdf",height=8,width=8)
plot(NA,ylim=c(6,8), xlim=c(1976,2005), ylab="Log GDP per capita", xlab="Time",axes=F)
lines(c(1976:2004),
      afripanel$lngdpmad[which(afripanel$Country=="Mauritania" &
                                 afripanel$year>="1976" 
                               & afripanel$year<="2004")],col="dark gray",lwd=2)
lines(c(1976:2004),
      afripanel$lngdpmad[which(afripanel$Country=="Niger" &
                                 afripanel$year>="1976" &
                                 afripanel$year<="2004")], col="dark gray", lwd=2)
lines(c(1976:2004), 
      afripanel$lngdpmad[which(afripanel$Country=="Chad" &
                                 afripanel$year>="1976" &
                                 afripanel$year<="2004")], col="dark gray", lwd=2)
lines(c(1976:2004),
      afripanel$lngdpmad[which(afripanel$Country=="Burundi" &
                                 afripanel$year>="1976" &
                                 afripanel$year<="2004")], col="black",  lty=2,lwd=3)
lines(c(1976:2004),
      afripanel$lngdpmad[which(afripanel$Country=="Ethiopia" &
                                 afripanel$year>="1976" &
                                 afripanel$year<="2004")], col="dark gray", lwd=2)
lines(c(1976:2004),
      afripanel$lngdpmad[which(afripanel$Country=="Guinea-Bissau" &
                                 afripanel$year>="1976" &
                                 afripanel$year<="2004")], col="dark gray", lwd=2)

lines(c(1976:2004), afripanel$lngdpmad[which(afripanel$Country=="Mali" &
                                               afripanel$year>="1976" &
                                               afripanel$year<="2004")], lwd=3)

abline(v=1991,lty=3, lwd=2,col="black") 
axis(1,at=c(1975,1980,1985,1990,1995,2000,2005))
axis(2,hadj=T)
legend("topleft",bty="n",legend=c("Mali","Burundi"),
       lty=c(1,2),lwd=c(3,3),cex=1,y.intersp=0.85)
graphics.off()


## Figure 1 -- Right panel
ind <- which(panel.reg$year>1975&
               !panel.reg$WBCode%in%c("BEN","GHA","KEN","MDG","MWI","SEN","TZA","ZAM"))
mali.mod <- (lm(lngdpmad ~ dem*poly(year,1)
                + lnpop
                + ki
                + openk
                + civwar
                + civwarend
                + wbank
                + imfadj
                + pwt_xrate
                + poly(year,1)*WBCode,
                panel.reg[ind,],y=T))

dat.t <- dat.c <- panel.reg[(rownames(panel.reg)%in%names(mali.mod$y)),]
dat.c$dem <- 0

treat <- predict(mali.mod,newdata=dat.t,type="response")[which(dat.t$WBCode=="MLI")]
control <- predict(mali.mod,newdata=dat.c,type="response")[which(dat.t$WBCode=="MLI")]

pdf(file="fig1-right.pdf",width=8,height=8)
plot(dat.t$year[which(dat.t$WBCode=="MLI")],dat.t$lngdpmad[which(dat.t$WBCode=="MLI")],
     type="l",ylim=c(6,8),ylab="Log GDP per capita",
     xlab="Time",lwd=3,bty="n",xlim=c(1976,2005))
lines(dat.t$year[which(dat.t$WBCode=="MLI")],control,lty=2,lwd=3)
legend("topleft",bty="n",legend=c("Mali","Fitted values under control"),
       lty=c(1,2),lwd=c(3,3),cex=1,y.intersp=0.85)
abline(v=1991,lty=3,lwd=2)
graphics.off()


## Figure 2
paths <- synth.wdk2.paths("MLI", "Mali", 1975, 1990, 1991, 2008, 6, 8)  

pdf(file="fig2.pdf",width=8,height=8)
plot(rownames(paths[[1]])[2:31],paths[[1]][2:31,1],
     type="l",ylim=c(6,8),ylab="Log GDP per capita",xlab="Time",lwd=3,bty="n")
lines(rownames(paths[[1]])[2:31],paths[[1]][2:31,2],lty=2,lwd=3)
legend("topleft",bty="n",legend=c("Mali","Synthetic counterfactual"),
       lty=c(1,2),lwd=c(3,3),cex=1,y.intersp=0.85)
abline(v=1991,lwd=2,lty=3)
graphics.off()

mali.weights <- round(paths[[2]]$solution.w,2)
country.names <- tapply(afripanel$Country,afripanel$wbcode2,unique)
rows <- unname(sapply(rownames(mali.weights),function(x){country.names[names(country.names)%in%x]}))
mali.weights <- data.frame(rows,mali.weights)
mali.weights <- mali.weights[order(mali.weights[,2],decreasing=T),]
mali.weights <- mali.weights[which(mali.weights[,2]!=0),]
colnames(mali.weights) <- c("Contribution","Weights")


## Figure 3
begin <- 1975
unitIDs <- unique(afripanel$WBCode)[unlist(tapply(
  afripanel$dem_ind,afripanel$WBCode,unique))==1][c(1:4,6,5,7:9)]
fullnames <- names(unlist(tapply(afripanel$dem_ind,afripanel$Country,unique)))[unlist(tapply(
  afripanel$dem_ind,afripanel$Country,unique))==1]
tr2s <- c(1989,1991,1992,1991,1993,1991,1993,1992,1991)
ends <- tr2s-1
final <- 2008
low <- 6
high <- 8
models.tr <- data.frame(unitIDs,rep(NA,length(unitIDs)))

for(i in 1:length(unitIDs)){
  models.tr[i,2] <- placebo.model.check(unitIDs[i],begin,ends[i],tr2s[i],final)
}
## Note: Error messages indicating non-conformable matrices may be printed when R
## fails to fit one of the two models encoded in the function "placebo.model.check"

results.tr <- data.frame(row.names = 1975:2008)
for(i in 1:length(fullnames)){
  #pdf(file=paste0(fullnames[i],".pdf"),width=8,height=8)
  if(models.tr[i,2]==1){results.tr[,i] <- synth.wdk1(unitIDs[i],fullnames[i],begin,ends[i],tr2s[i],final,low,high)}
  else{results.tr[,i] <- synth.wdk2(unitIDs[i],fullnames[i],begin,ends[i],tr2s[i],final,low,high)}
  graphics.off()
}
colnames(results.tr) <- unitIDs

alphas <- c()
negs <- c()
for(i in 1:ncol(results.tr)){
  alphas[i] <- (mean(results.tr[which(rownames(results.tr)==(tr2s[i]+1)):nrow(results.tr),i]))
  negs[i] <- mean(sapply(results.tr[which(rownames(results.tr)==(tr2s[i]+1)):nrow(results.tr),i],
                         function(x){ifelse(x<0,1,0)}))
}
names(alphas) <- names(negs) <- colnames(results.tr)

mspe.prior <- c()
for(i in 1:ncol(results.tr)){
  mspe.prior[i] <- (mean(results.tr[1:which(rownames(results.tr)==(tr2s[i])),i]^2))
}

pdf(file="fig3-benin.pdf")
synth.wdk1("BEN", "Benin", 1975, 1988, 1989,2008,6,8)
text(labels=paste("MSPE","=",round(mspe.prior[1],4),sep=" "),x=mean(c(1975,tr2s[1])),y=7.75)
graphics.off()

pdf(file="fig3-ghana.pdf")
synth.wdk2("GHA", "Ghana", 1975, 1990, 1991, 2008, 6, 8)
text(labels=paste("MSPE","=",round(mspe.prior[2],4),sep=" "),x=mean(c(1975,tr2s[2])),y=7.75)
graphics.off()

pdf(file="fig3-kenya.pdf")
kenya <- synth.wdk2("KEN", "Kenya",1975, 1991, 1992, 2008,6, 8) 
text(labels=paste("MSPE","=",round(mspe.prior[3],4),sep=" "),x=mean(c(1975,tr2s[3])),y=7.75)
graphics.off()

pdf(file="fig3-madagascar.pdf")
synth.wdk2("MDG", "Madagascar", 1975, 1991, 1992, 2008, 6, 8) 
text(labels=paste("MSPE","=",round(mspe.prior[4],4),sep=" "),x=mean(c(1975,tr2s[4])),y=7.75)
graphics.off()

pdf(file="fig3-malawi.pdf")
synth.wdk1("MWI", "Malawi", 1975, 1992,1993,2008, 6, 8)
text(labels=paste("MSPE","=",round(mspe.prior[5],4),sep=" "),x=mean(c(1975,tr2s[5])),y=7.75)
graphics.off()

pdf(file="fig3-mali.pdf")
synth.wdk2("MLI", "Mali", 1975, 1990, 1991, 2008, 6, 8)
text(labels=paste("MSPE","=",round(mspe.prior[6],4),sep=" "),x=mean(c(1975,tr2s[6])),y=7.75)
graphics.off()

pdf(file="fig3-senegal.pdf")
synth.wdk2("SEN", "Senegal",1975, 1992, 1993, 2008, 6, 8) 
text(labels=paste("MSPE","=",round(mspe.prior[7],4),sep=" "),x=mean(c(1975,tr2s[7])),y=7.75)
graphics.off()

pdf(file="fig3-tanzania.pdf")
synth.wdk1("TZA", "Tanzania", 1975, 1991,1992,2008,6, 8) 
text(labels=paste("MSPE","=",round(mspe.prior[8],4),sep=" "),x=mean(c(1975,tr2s[8])),y=7.75)
graphics.off()

pdf(file="fig3-zambia.pdf")
zambia <- synth.wdk1("ZMB", "Zambia", 1975, 1990, 1991, 2007, 6, 8)
text(labels=paste("MSPE","=",round(mspe.prior[9],4),sep=" "),x=mean(c(1975,tr2s[9])),y=7.75)
graphics.off()


##Figure 4
results.tr1 <- results.tr[,order(apply(results.tr[which(rownames(results.tr)=="1991"):
                                                    nrow(results.tr),],2,sum))]

pdf(file="fig4.pdf")
par(mai=c(1.2,0.9,0.5,0.5),oma=c(0,0,0,0))
plot(NA, xlim=c(-0.4,0.6), ylim=c(0,9.5), axes=F, ylab="",
     xlab="Distributions of annual logged post-treatment differences in per capita GDP (USD)")
for(i in 1:ncol(results.tr)){
  boxplot(results.tr1[which(rownames(results.tr1)==(tr2s[i]+1)):nrow(results.tr1),i],
          horizontal=T,add=T,axes=F,at=i,outline=F)
}
abline(v=0,lty=2,lwd=2)
text(y = seq(1, 9, by=1), x=par("usr")[1], par("usr")[3] - 0.2,
     labels = (c("Zambia","Madagascar","Tanzania","Malawi","Kenya","Benin","Ghana",
                 "Senegal","Mali")),
     cex=1,srt = 45, pos = 2, xpd = TRUE)

axis(1)
graphics.off()


## Table 1
set.seed(1) ## Set randomization seed
## Pooled OLS
res0 <- c()
N <- 1000
for(i in 1:N){
  blocks <- sample(unique(panel.reg$wbcode2),length(unique(panel.reg$wbcode2)),replace=T)
  s <- as.vector(sapply(blocks,function(x){which(panel.reg$wbcode2%in%x)}))
  mod <- lm(lngdpmad ~ dem
            + lnpop
            + ki
            + openk
            + civwar
            + civwarend
            + wbank
            + imfadj
            + eximdiff
            + pwt_xrate
            + year
            , data=panel.reg[s,]
            , y=T
  )
  res0[i] <- coefficients(mod)[2]
}

## Hausman Test -- Country-specific intercepts
fe.mod <- lm(lngdpmad ~ dem
             + lnpop
             + ki
             + openk
             + civwar
             + civwarend
             + wbank
             + imfadj
             + eximdiff
             + pwt_xrate
             + year
             + factor(wbcode2)
             , data=panel.reg
             , y=T
)

re.mod <- lmer(lngdpmad ~ dem
               + lnpop
               + ki
               + openk
               + civwar
               + civwarend
               + wbank
               + imfadj
               + eximdiff
               + pwt_xrate
               + year
               + (1|wbcode2)
               , data=panel.reg
)

betadiff <- coefficients(fe.mod)[2:12] - fixef(re.mod)[2:12]
vardiff <- vcov(fe.mod)[2:12,2:12] - vcov(re.mod)[2:12,2:12]
H <- t(betadiff) %*% solve(vardiff) %*% betadiff
ifelse(pchisq(as.numeric(H), df=11, lower.tail=F)<=0.05,("Reject the null (use FE)"),
       ("Fail to reject the null (use RE)"))

## Country-random intercepts
res1 <- c()
N <- 1000
for(i in 1:N){
  blocks <- sample(unique(panel.reg$wbcode2),length(unique(panel.reg$wbcode2)),replace=T)
  s <- as.vector(sapply(blocks,function(x){which(panel.reg$wbcode2%in%x)}))
  mod <- lmer(lngdpmad ~ dem
              + lnpop
              + ki
              + openk
              + civwar
              + civwarend
              + wbank
              + imfadj
              + eximdiff
              + pwt_xrate
              + year
              + (1|wbcode2)
              , data=panel.reg[s,]
              #, y=T
  )
  res1[i] <- fixef(mod)[2]
}

## Hausman Test -- Country-varying slopes and time-trends
fe.mod <- lm(lngdpmad ~ dem
             + lnpop
             + ki
             + openk
             + civwar
             + civwarend
             + wbank
             + imfadj
             + eximdiff
             + pwt_xrate
             + factor(year)
             + factor(wbcode2)
             , data=panel.reg
             , y=T
)

re.mod <- lmer(lngdpmad ~ dem
               + lnpop
               + ki
               + openk
               + civwar
               + civwarend
               + wbank
               + imfadj
               + eximdiff
               + pwt_xrate
               + (1|year)
               + (1|wbcode2)
               , data=panel.reg
)

betadiff <- coefficients(fe.mod)[2:11] - fixef(re.mod)[2:11]
vardiff <- vcov(fe.mod)[2:11,2:11] - vcov(re.mod)[2:11,2:11]
H <- t(betadiff) %*% solve(vardiff) %*% betadiff
ifelse(pchisq(as.numeric(H), df=11, lower.tail=F)<=0.05,("Reject the null (use FE)"),
       ("Fail to reject the null (use RE)"))

## Country-varying intercepts and time trends
res2 <- c()
N <- 1000
for(i in 1:N){
  blocks <- sample(unique(panel.reg$wbcode2),length(unique(panel.reg$wbcode2)),replace=T)
  s <- as.vector(sapply(blocks,function(x){which(panel.reg$wbcode2%in%x)}))
  mod <- lm(lngdpmad ~ dem
            + lnpop
            + ki
            + openk
            + civwar
            + civwarend
            + wbank
            + imfadj
            + eximdiff
            + pwt_xrate
            + factor(year)
            + factor(wbcode2)
            , data=panel.reg[s,]
            , y=T
  )
  res2[i] <- coefficients(mod)[2]
}

## Country-varying slopes and intercepts
res3 <- c()
N <- 1000
for(i in 1:N){
  blocks <- sample(unique(panel.reg$wbcode2),length(unique(panel.reg$wbcode2)),replace=T)
  s <- as.vector(sapply(blocks,function(x){which(panel.reg$wbcode2%in%x)}))
  mod <- lm(lngdpmad ~ dem*factor(wbcode2)
            + lnpop*factor(wbcode2)
            + ki*factor(wbcode2)
            + openk*factor(wbcode2)
            + civwar*factor(wbcode2)
            + civwarend*factor(wbcode2)
            + wbank*factor(wbcode2)
            + imfadj*factor(wbcode2)
            + eximdiff*factor(wbcode2)
            + pwt_xrate*factor(wbcode2)
            + year*factor(wbcode2)
            , data=panel.reg[s,]
            , y=T
  )
  res3[i] <- coefficients(mod)[2]
}

panel.estimates <- cbind(c(mean(res0),mean(res1),mean(res2),mean(res3)),
      c(ecdf(res0)(0),ecdf(res1)(0),ecdf(res2)(0),ecdf(res3)(0)))
rownames(panel.estimates) <- c("Pooled OLS","Country random effects","Country and year fixed effects",
                               "Country-varying slopes and intercepts")
colnames(panel.estimates) <- c("Estimates","P-values")


## Table 2
moderators <- matrix(data=c(cor.test(a$inst_sys,a$effect)$estimate,
              cor.test(a$inst_sys,a$effect)$p.value,
              -cor.test(a$unicameral,a$effect)$estimate,
              cor.test(a$unicameral,a$effect)$p.value,
              cor.test(ifelse(a$pr==0,1,0),a$effect)$estimate,
              cor.test(ifelse(a$pr==0,1,0),a$effect)$p.value,
              cor.test(ifelse(a$protest==0,1,0),a$effect)$estimate,
              cor.test(ifelse(a$protest==0,1,0),a$effect)$p.value,
              cor.test(a$gini,a$effect)$estimate,
              cor.test(a$gini,a$effect)$p.value,
              cor.test(a$econ_lib_prior,a$effect)$estimate,
              cor.test(a$econ_lib_prior,a$effect)$p.value),
       nrow=6,ncol=2,byrow=T)
rownames(moderators) <- c("Presidentialism","Bicameralism","Plurality voting",
                               "Transition under foreign pressure","Income inequality",
                               "Economic liberalization before democratization")
colnames(moderators) <- c("Estimates","P-values")


## Figure 5
placebo <- c(unique(afripanel$WBCode[afripanel$cont_dem_ind==1]))

start <- 1975; end <- 1991; tr2 <- 1992; final <- 2005

models <- data.frame(placebo,rep(NA,length(placebo)))
for(i in 1:length(placebo)){
  if(placebo[i]=="ETH"|placebo[i]=="SDN"){next}    
  models[i,2] <- placebo.model.check(placebo[i],start,end,tr2,final)
}
models[is.na(models[,2]),2] <- 1    ## Data limitations model 1 must be used for Ethiopia and Sudan

start <- 1975
gaps.plac <- matrix(NA, nrow=(final-start)+1, ncol=length(placebo))
for(i in 1:length(placebo)){
  if(models[i,2]==1){spe <- synth.placebo1(placebo[i],start,end,tr2,final)}
  else{spe <- synth.placebo2(placebo[i],start,end,tr2,final)}
  gaps.plac[,i] <- spe
}

rownames(gaps.plac) <- start:final
colnames(gaps.plac) <- placebo

spe.plac <- gaps.plac^2

mspe.plac <- matrix(NA,1,length(placebo))
for (i in 1:length(placebo)){
  mspe.plac[,i] <- mean(spe.plac[which(rownames(spe.plac)==start):which(rownames(spe.plac)==end),i])
}

gaps.plac <- sort(apply(gaps.plac[(end-start):(final-start),which(mspe.plac<0.01)],2,sum))
gaps.treat <- sort(apply(results.tr[(end-start):(final-start),],2,sum))

graph <- melt(list(gaps.treat,gaps.plac))
graph[,3] <- 1
graph[c(10,11,12,18),3] <- 2     ## "Placebo - Civil War" cases include Angola, Rwanda, Ethiopia,
                                 ## and Sudan (Themner and Wallensteen, 2013)
graph[c(13:17,19:22),3] <- 3

colnames(graph) <- c("value","treated","alpha")
rownames(graph) <- c(names(gaps.treat),names(gaps.plac))

labs <- c("Zambia (-)","Madagascar (-)","Tanzania (-)","Kenya (nil.)",
          "Malawi (nil.)","Senegal (+)","Benin (+)","Ghana (+)","Mali (+)")

xpos <- sapply(graph[which(graph[,2]==1),1],function(y){if(y<0)
               {(round_any(y,accuracy=0.3,f=floor))}
               else {(round_any(y,accuracy=0.3,f=round))}})
xpos <- xpos + ifelse(xpos[1:9]<=0,0.15,-0.15)

ypos <- c(ifelse(diff(xpos)==0,1.5,0.5),0.5)

pdf(file="fig5.pdf",width=6.78,height=3.5)
ggplot(graph, aes(value, fill = as.factor(alpha))) + 
  geom_histogram(binwidth=0.3) + 
  scale_fill_manual(limits=c(1,3,2),breaks=c(1,3,2),values=c("black","gray50","gray80"),
                    guide=guide_legend(title=""),
                    labels=c("Treated","Placebo","Placebo - Civil War")) + 
  theme_bw() + 
  theme(axis.title=element_text(size=8),axis.text=element_text(size=8),
        legend.text=element_text(size=8),legend.key.size=unit(0.5, "cm")) + 
  annotate("text",x=xpos, y=ypos, label=labs, angle=90, color="white", size=2) + 
  xlab("Cumulative Post-Treatment Difference in log(GDP per Capita)") + 
  ylab("Frequency")
graphics.off()

## Save the console output as "growth-dem-wdk-tables.txt" to working directory
sink(file="growth-dem-wdk-tables.txt",type="output")

## Figure 2 table
print(mali.weights,row.names=F)

## Table 1
panel.estimates

## Table 2
moderators

sink()